home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SVDFIT.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  75 lines
  1. PROGRAM d14r4(input,output);
  2. (* driver for routine SVDFIT *)
  3. (* polynomial fit *)
  4. CONST
  5.    npt=100;
  6.    spread=0.02;
  7.    npol=5;
  8.    mp=npt;
  9.    np=npol;
  10.    ncvm=npol;
  11. TYPE
  12.    glndata = ARRAY [1..npt] OF real;
  13.    glmma = ARRAY [1..npol] OF real;
  14.    glnpbynp = ARRAY [1..np,1..np] OF real;
  15.    glmpbynp = ARRAY [1..mp,1..np] OF real;
  16.    glnparray = glmma;
  17.    glmparray = glndata;
  18.    (* The previous 2 declarations are required because of the absence *)
  19.    (* of conformant arrays in most Pascal implementations. *)
  20.    glcvm = ARRAY [1..ncvm,1..ncvm] OF real;
  21. VAR
  22.    glinext,glinextp : integer;
  23.    glma : ARRAY [1..55] OF real;
  24.    gliset : integer;
  25.    glgset : real;
  26.    chisq : real;
  27.    i,idum : integer;
  28.    x,y,sig : glndata;
  29.    a : glmma;
  30.    cvm : glcvm;
  31.    u : glmpbynp;
  32.    v : glnpbynp;
  33.    w : glnparray;
  34.  
  35. (*$I MODFILE.PAS *)
  36. (*$I RAN3.PAS *)
  37.  
  38. (*$I GASDEV.PAS *)
  39.  
  40. (*$I SVDCMP.PAS *)
  41.  
  42. (*$I SVBKSB.PAS *)
  43.  
  44. (*$I SVDVAR.PAS *)
  45.  
  46. PROCEDURE func(x: real; VAR p: glmma; mma: integer);
  47. (* This is essentially FPOLY renamed. *)
  48. VAR
  49.    j: integer;
  50. BEGIN
  51.    p[1] := 1.0;
  52.    FOR j := 2 to mma DO p[j] := p[j-1]*x
  53. END;
  54.  
  55. (*$I SVDFIT.PAS *)
  56.  
  57. BEGIN
  58.    gliset := 0;
  59.    idum := -911;
  60.    FOR i := 1 to npt DO BEGIN
  61.       x[i] := 0.02*i;
  62.       y[i] := 1.0+x[i]*(2.0+x[i]*(3.0+x[i]*(4.0+x[i]*5.0)));
  63.       y[i] := y[i]*(1.0+spread*gasdev(idum));
  64.       sig[i] := y[i]*spread
  65.    END;
  66.    svdfit(x,y,sig,npt,a,npol,u,v,w,mp,np,chisq);
  67.    svdvar(v,npol,np,w,cvm,npol);
  68.    writeln;
  69.    writeln('Polynomial fit:');
  70.    FOR i := 1 to npol DO BEGIN
  71.       writeln(a[i]:12:6,'  +-',sqrt(cvm[i,i]):10:6)
  72.    END;
  73.    writeln('Chi-squared',chisq:12:6);
  74. END.
  75.